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A generalized form of the parametric 
spline methods of degree (2k + 1) for 
solving a variety of two-point boundary 
value problems 


Z. Sarvari 


Abstract 


In this paper, a high order accuracy method is developed for finding the 
approximate solution of two-point boundary value problems. The present 
approach is based on a special algorithm, taken from Pascal’s triangle, for 
obtaining a generalized form of the parametric splines of degree (2k + 1), 
k =1,2,..., which has a lower computational cost and gives the better ap- 
proximation. Some appropriate band matrices are used to obtain a matrix 
form for this algorithm. 

The approximate solution converges to the exact solution of order 
O(h4*), where k is a quantity related to the degree of parametric splines 
and the number of matrix bands that are applied in this paper. Some 
examples are given to illustrate the applicability of the method, and we 
compare the computed results with other existing known methods. It is 
observed that our approach produced better results. 
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1 Introduction 


Spline function is a piecewise polynomial satisfying certain conditions of the 
continuity of the function and its derivatives. In other words, a spline function 
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S(x) of degree d is defined in a region [a,b] such that there exists a mesh 
A = {a= 2% < 4 < fg <0) < Un < Up = Db} with h; = x; — x;_1 for 


i=1,2,...,n. This function satisfies the following conditions: 
(i) In each subinterval [v;, 7:41], 7 = 0,1,...,2—1, S(x) is a polynomial 
of degree d. 


(ii) S(x) and its first (d — 1) derivatives are continuous on [a, }]. 

The spline’s theory and application were thoroughly discussed by Ahlberg 
Nilson, and Walsh [1] and Greville [9]. So far, different types of spline meth- 
ods, such as approximating, interpolating, and curve fitting functions, have 
been developed and used to solve a wide variety of differential equations; see, 
for example, [2, 3, 14, 15, 17, 20, 22, 26, 25, 24, 31] and references therein. 
One type of splines, considered in this paper, is the parametric splines de- 
veloped to address some shortcomings of ordinary spline methods. These 
splines, depending on a parameter 7 > 0, are defined through the solution 
of a differential equation in each subinterval. The arbitrary constants are 
chosen to satisfy certain smoothness conditions at the joints. These splines 
reduce to polynomial splines (ordinary splines) as 7 > 0. The exact form 
depends upon the manner in which the parameter is introduced. Therefore, 
different types of parametric splines with distinct convergence orders can be 
generated. Although, these methods obtained significant results, due to the 
lengthy calculations, no attempt was made to extend parametric splines of 
higher degrees. Note that using the word “degree”, in this paper, for para- 
metric splines is only for numbering and ordering them. It does not have the 
common meaning that is used for polynomials. 

For the first time, this paper presents a general form of parametric splines 
with the degree (2k +1), k =1,2,..., which has a lower computational cost 
and a higher-order of convergence than the usual methods using parametric 
splines. Before going into details about the method, it seems necessary to 
review some of the fundamental properties and definitions of these parametric 
splines in the following subsection. 


1.1 Parametric spline methods with the degree (2k + 1) 


For simplicity, it is assumed that the subintervals are of equal length, so 
h = h; = hii. The interval [a,b] is divided into n equal subintervals using 
knots a; and the partition A = {a = 20, 21,...,@n = b}, where x; = zo + th, 
h = °=* and n is a positive integer. The parametric spline function S(2), 
with the degree (2k +1), k = 1,2,..., is obtained in the subinterval [x;_1, xj] 
by solving the following differential equation and determining the constants 
of integration: 


TL — Xj-1 


ger) (x) ai 72 S(2k-2) (7) = (S*) (x;) fe 72g (2k—-2) (x;))( ) 
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This function of class C?*[a,b] depends on a parameter T and reduces 
to an ordinary spline function with the degree (2k +1), as tr + 0. The 
continuity of its derivatives at the grid points, that is, Ss (2%) - 3) (xi), 
vy =1,3,...,2k-—1, yields spline relations. Note that 5; is the spline function 
in the subintervals [x;, 7:41]. Using algebraic manipulation on these relations, 
a differential relation, called “consistency relation”, is obtained in terms of u 
and its derivatives at knots. In the parametric spline methods, the approx- 
imate solution of a given boundary value problem (BVP) is determined by 
solving the system defined by this consistency relation. 

Now, to further explain how parametric splines are used to solve equa- 
tions, we consider a simple second order BVP as follows: 


u(x) = f(a) + g(a)ju(z), aw € [a, B], 
u(a) = A, (1) 
u(b) = 7, 


where \ and 7¥ are finite real constants and the functions f(x) and g(x) are 
continuous on [a,b]. Such problems arise in the theory that describes the 
deflection of plates and a number of other scientific applications [10]. 

The consistency relation associated with (1) in spline methods, is in terms 
of wu; and ul’. Note that u; = u(#;) and ul’ = u"(x,). A system of linear alge- 
braic equations is generated by substituting discretized (1) in the mentioned 
consistency relation. Finally, by solving this system, the approximate so- 
lution of (1) is obtained. One can observe that, for k > 1, the number of 
equations in this system is less than the number of unknowns; see, for exam- 
ple, [2, 3, 4, 8, 7, 11, 13, 14, 16, 18, 20, 19, 21, 26, 25, 24, 30, 31] To obtain 
the unique solution of the system, more equations, called “end conditions or 
boundary formulas”, are needed. 

When & is a large number, two problems are encountered. First, the 
number of additional equations that must be defined to complete the afore- 
mentioned system increases. Second, as the value of k increases, so does the 
number of relations resulting from the derivative continuities of splines, and 
consequently, the combination of them becomes more difficult. 

In this paper, we present a method that allows us to obtain a general form 
for the consistency relations of parametric spline methods of degree (2k + 1) 
without going through a lengthy and complex calculation process. For large 
values of k, it does not face the drawbacks mentioned above, because, in the 
proposed method, we do not need to obtain the spline function and use its 
continuity properties and derivatives directly. This means that we do not 
solve any differential equation. In fact, by providing a general pattern, both 
the consistency relation and the required additional equations are obtained 
without a need to generate many spline relations. Moreover, we transform 
the desired algorithm into a matrix form by defining proper band matrices, 
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which gives us more insight into the method and facilitates convergence study. 
Furthermore, the convergence order of splines is improved by this general 
formulation. 


We apply our method to (1), which is in terms of u and wu”. However, the 
method can be applied to more complex models of (1), such as the nonlinear 
form or the system of these equations. It will be demonstrated in the nu- 
merical results section. Our method along with Newton—Raphson method is 
used to solve Bratu’s problem (Example (2)) as a nonlinear equation. Also, 
our method is applied to solve problems such as Perturbed (Example (3)) 
and two-dimensional problems in the calculus of variations (Example (4)). 


It should be noted that, if an equation other than (1) is considered, then 
the execution process of the method, such as generating a consistency rela- 
tion and additional equations, will be changed. Because they are produced 
and defined according to the type of equation. In addition, while this paper 
focuses on parametric splines with the degree (2k + 1), the extending of our 
method can be probed for other types of splines as well, that is, nonpoly- 
nomials, ordinary splines, and parametric splines with the degree of (2k), 
k = 2,3,4,.... 

The outline of this paper is as follows. In section 2, a comprehensive 
description of the method is given. To demonstrate the efficiency and su- 
periority of the presented method, we solve examples of linear, nonlinear, 
perturb, and system of two-point BVPs and compare the obtained results 
with the other quoted methods in section 3. Finally, some important con- 
cluding remarks are given in section 4. 


2 Derivation of the method 


In this section, we describe our method in detail. As mentioned previously, 
in the common form of the spline method, we need to generate a consistency 
relation proportional to the type of BVP that we are going to solve numeri- 
cally. This can be time-consuming and even complicated due to the numerous 
calculations required, such as developing the spline function criterion, deter- 
mining its coefficients, and computing successive derivatives. We provide a 
generalized form for the consistency relation of all parametric spline methods 
of degree (2k+1), while solving (1), without the need for long calculations. To 
produce this relation, we find a specific pattern and then convert it to a matrix 
form by using only the properties of band matrices, especially, the following 
widely used a matrix C, which is (n — 1) x (n — 1)-dimensional and evident 
in the majority of spline-based papers (see [8, 7, 13, 14, 18, 20, 19, 21, 31]): 


2, t=), 
(C)ey = 4-1, | =| =1, (2) 
0, otherwise. 
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To shed light on the above-mentioned issues, we begin subsection 2.1 
by evaluating some samples of the consistency relations associated with the 
spline methods previously used by researchers and then identifying a general 
pattern for our desired consistency relation. In subsection 2.2, we define its 
matrix form. Subsections 2.3 and 2.4 are also dedicated to solving (1) and 
developing boundary formulas according to the contents of the previous two 
subsections. 


2.1 The consistency relation 


There are two types of coefficients in the consistency relations of spline meth- 
ods: the coefficients of wu and its derivatives. By studying the spline papers 
(see the references on page 3), we find that just the coefficients of u follow a 
certain pattern. The coefficients of u are of two kinds: known and unknown. 
The first type of coefficients exists in the consistency relations of splines with 
the degree (24 + 1) in solving the BVPs of order (2k), while the second ones 
can be seen in the consistency relations of the same splines in solving the 
BVPs of order 2,4,6,...,2k — 2, for k = 2,3,4,... (ie., for k = 2, we have 
spline with the degree 5 in solving BVP of order 2, for k = 3, we have spline 
with the degree 7 in solving BVP of order 2 and 4, for k = 4, we have spline 
with the degree 9 in solving BVP of order 2,4 and 6, etc.) We first establish 
a pattern for known coefficients and then extend this to unknown ones. 

In the following, we highlight them in the sample format. For convenience, 
we consider the coefficients of u to be on the left side of the consistency 
relation and assume that the coefficients of derivatives of u be on the right 
side. 

Sample 1 (k=1): The left side of the consistency relation of spline 
method with the degree 3 in solving a BVP of order 2: 

For 1=1,2,...,n—1: 


luj_1 = 2u; + Luj+i Sree, 


One can see this relation in [15, 29]. 

Sample 2 (k=2): The left side of the consistency relation of spline 
method with the degree 5 in solving a BVP of order 4: 

For 7 = 2,3,...,n—2, 


Luj_2 = Auj_1 + 6u; = Auisi + Lluj+e ae ER 


One can see this relation in [18]. 

Sample 3 (k=3): The left side of the consistency relation of spline 
method with the degree 7 in solving a BVP of order 6: 

For 7 = 3,4,...,n —3, 
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Luj;_3 — 6uj_2 + 15u,;_1 — 20u; + 15uj;41 — 6uj4e + luja3 =:--. 


One can see this relation in [3, 11, 26]. 

Sample 4 (k=4): The left side of the consistency relation of spline 
method with the degree 9 in solving a BVP of order 8: 

For 1 = 4,5,...,n—A4, 


luj_4 = 8uj_3 + 28uj;_2 as 56uj;_1 + 70u; 
— 56uj41 + 28ui42 — 8ui43 + Luiza =-::. 

One can see this relation in references [2, 19]. 

Sample 5 (k=5): The left side of the consistency relation of spline 
method with the degree 11 in solving a BVP of order 10: 

For 7 = 5,6,...,n—5, 

luj_s a 10u;_4 =) 45uj;_3 _ 120u;_2 + 210u;_1 = 252u,; 

Tr 210uUi+1 — 120u;+2 + 45ui43 — 10uj+4 + Lluiis Se 


One can see this relation in [20, 24]. 

Sample 6 (k=6): The left side of the consistency relation of spline 
method with the degree 13 in solving a BVP of order 12: 

For 1 = 6,7,...,n—6, 


lu;_¢ — 12u;_5 + 66u;_4 — 220u;_3 + 495u;_2 — 792u;_; + 924u; 
— T92ui41 + 495uUji42 — 220Ui43 + 66uj+4 — 12Uj45 +1luiy6 =--:- 


One can see this relation in [25]. 

By considering the above coefficients, we find that, regardless of their 
sign, they are the same as the binomial coefficients or the entries in the rows 
of Pascal’s triangle: 


121 
1 3.3 1 
146441 
1 5 10 10 5 1 
1 6 15 20 15 6 1 
1 7 21 35 35 21 7 1 
1 8 28 56 70 56 28 8 1 
1 9 36 84 126 126 84 36 9 1 
1 10 45 120 210 252 210 120 45 10 1 
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1 11 55) 165 330 462 462 330 165 55 «11 #1 
1 12 66 220 495 792 924 792 495 220 66 12 1 


The correlation between the above-known coefficients of u and Pascal’s 
triangle motivates us to find a similar correlation for the unknown ones that 
we will deal with in this paper. 

After studying the references such as [8, 14, 15, 16, 21, 22, 29, 30, 31] (Pre- 
vious studies have only investigated 3rd-, 5th-, 7th-, and 9th-degree spline 
methods. Higher degree splines have not been used yet), we find out to con- 
sider the initial form for the consistency relations of spline methods with the 
degree (2k + 1), k = 2,3,4,... in solving a BVP of order two as follows: 


* (Wik + Ute) + * (Uinta + UWipe—1) $e + (Usa + Ui41) 
+ ui = —h? (Boul + Bi(uj_y + fer) +--+ Balu, + Ue) » 


where {;’s are the coefficients which will be determined numerically during 
the process of the method. We have displayed the vacancy of the unknown 
coefficients of wu with *. We intend to find a pattern for them, inspired 
by Pascal’s triangle. For this purpose, we first consider the parameters as 
Qo, @1,2,---,Az—1, for each k, and then we implement Pascal’s algorithm for 
them. It should be mentioned that the numerical value of these coefficients 
for each k is independent of the values for other k, so it is preferable to write 
6;’s and a,;’s with the exponent (k) as eo and al”, However, to reduce 
the complexity of the text, the exponent (k) could be removed from the 
coefficients without disturbing the whole. Indeed, the number of coefficients, 
which is indicated by an index in them, is affected by k. 


Hence, we have the following Pascal’s algorithm for ag, a1, @2,...,Qg—1: 
Ak—-1;,%k—2,--+5 2, A1,A0,%1,02,---, Ak—-2,%k-1 
O~—1, &h—1 + Ap—2,---; a, +a9,09 + 4,---, A~—-2 + Ap-1,%k-1 
Op—1, 204-1 + Ap—2,---, a2 + 2a1 + a0, 2(a1 + a0), a2 + 201 + a0,.--, 2p, -1 + Ap—2,%K-1 


According to the consistency relations of the mentioned references, the 
third row of the above triangle, regardless of the signs, is the same as the 
coefficients of u in the consistency relation of the spline method of degree 
(2k + 1) for solving a BVP of second order as (1). We will demonstrate 
that the next rows of this triangle are the coefficients of u in the consistency 
relations of the spline methods of degree (2k + 1) for solving BVPs of higher 
orders (bigger than 2) in future research. 


Thus, the following relation with the sign (—1)?*? for each phrase apuitg, 
can be defined as the desired consistency relation for k = 2,3,4,...: 
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—Op—1(Ui-k + With) + (2aK—1 — Ak—2)(Ui-k41 + Uitk—1) 
+(—ap—-1 + 2aK—2 — An—3) (Wik 42 + UstR—2) + +++ + (200 — 201) Ui 


k 
=—h? | Boul +> Bi(ui_;+uiys) |, t=h k+1,....n-k. (3) 
j=l 


In the following, to verify the correctness of (3), we compare it to the 
consistency relations developed in related papers. A quick review shows that 
although the appearance of the coefficients in the consistency relations of 
available references slightly differs from what we propose, they are identical 
in content. In fact, the other authors have defined these coefficients in terms 
of parameter T: 


Equation (3) for k = 2, is the consistency relation of parametric 
spline with degree 5 (quintic spline): 


— a1 (Uj—2 + Ui42) + (2Za1 — ag) (us—1 + Wi41) + (— 201 + 200) Ui 
= —h? [B2(ut_» a Ups) a5 Bi (uj_y + Uy41) 7 Bout | ? i= 2, 3, ao 2. 
One can compare it to the consistency relations in [16, 21, 30, 31]. 


Equation (3) for k = 3, is the consistency relation of parametric 
spline with degree 7 (septic spline): 
—a2(uj_3 + Uit3) + (2a = Q1)(uj_2 + Ui+2) 
+(—ao + 2a; — Q2)(Ui—1 + Ui41) + (—2a, + 200) ui 
=—h? [83 (uj_3 + Uy,3) + Bo(uj_> + Uj) + By (ujy_y + Uy 41) + Bou;], 
1=3,4,...,n—3. 


One can compare it to the consistency relations in [14]. 


Equation (3) for k = 4, is the consistency relation of parametric 
spline with degree 9 (nonic spline): 


—a3z(uja + Ui44) + (203 — a2)(ui-3 + Ui43) + (—a3 + 2a — a1)(ui—2 + Ui42) 
+(—ag + 2a1 — ao) (ui-1 + Ui41) + (—201 + 209) Uj 
= —h?[Ba(uj_4 + wpa) + Bs(uy_g + wig) + Bo(uj_o + wip), 
+B, (ui, + ui) + Bou;] 4=4,5,...,n—4. 


One can compare it to the consistency relations in [8, 7]. 
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2.2 The matrix form 


To demonstrate the accuracy of the above cases, namely, the validity of our 
claim about the correlation between the consistency relation of the parametric 
spline of degree (2k + 1) and Pascal’s algorithm, we need to obtain a matrix 
form for (3), dependent on k. For this purpose, we use a matrix C and the 
following band matrices, which are (n — 1) x (n — 1)-dimensional and play a 
major role in our method: 


Qo, t= J, 
1, |i — jl = 1, 
(Ajig = 4: (4) 
Qk-1; ji-—j)/ =k-1, 
0, otherwise, 
Bo, t=J, 
Pi, |i _ JI = 1, 
(Big =4' 


Br, |i _ j| = k, 
0, otherwise. 


Since the coefficients of the mentioned samples of consistency relations 
can be observed in the rows of the consecutive multiplication of matrix C by 
itself, it is expected that the coefficients of (3) is obtained from multiplying 
matrix A by the matrix C. Therefore, for each k = 1,2,..., the general 
matrix form can be defined for (3) as follows: 


(AC) paieV +h? (B)egi~V” =0, 4 =k, k+1,...,.n—&. (5) 


where 4,41, denotes (& + 1)th row of the above matrices. Proportional to k, 
we also define the column vectors V and V” as follows: 


V, = d Mir kt+G-1) 1<j <2k+1, 
7 ho; MAk+1<j<n-l, 
and 
via) ey 1 STs et) 
a 0, 2ke+1<j<n-l, 


where V; and V;’ are the jth element of vectors V and V", respectively. 
Now, according to (5), we can obtain a set of parametric spline methods 

by changing &, without need to know the nature of the elements of matrices A 

and B. The numerical values of these elements are determined by expanding 
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(5) in Taylor’s series around «;. For this purpose, we first rewrite (5) in the 


following form: 


(AC) n41,1(Ui-k + With) + (AC)n41,2(ui-k41 + Uite—1) +... 


I 


+++++ (B)psipyiu;) 
= ee Se ee ee 


Note that in the above relation, we used the equalities 


+(AC) 41,4410 + A? ((B) e411 (uy + fee) + (B)egt2(uf naa + Uf e_1) 


(AC)k41,5 = (AC)e41,2¢+(2-7) and (B)esiy = (B)e+ire+e-j, 1S ji < ke. 
These are directly obtained from the definition of band matrices A, B and 


C. 


Now the local truncation error T;, corresponding to Taylor’s series of (5) 


can be obtained as 


(—kh)? ” 
2! 


kh)? 
T; =(AC)e41,1(ui — khuj + (kh) 


2! 
= 2 
(RADA yy 


+ (AC) 41,2(us + (—k + 1)hu; 4 


k —1)h)? 
+ uz + (k- 1)hui 4 (( = uy f... Hee + (AC) p41 pi 


—kh)? 
+ h?((B)kqra(uy — khu® +! mh ul) 4... ul! + khu® 


ue) 4.) + (B)eg2(uy + (-k + hu? + 


((k = 1h)? 
2! 


2! 


ull + (k—1)hul? + Tey © ane 


+ (B)k-+146+14j )- 
On simplifying, we get 


T; = (2(AC) e411 + 2(AC) pgtje +++ + 2(AC)egik + (AC) e412, 041) Ui 


+2(B)ipi1 + 2(B)epie t++> + 2(B)etie + (Bye tiepih7uy 


Ul ees tu + khu, + ul +... 


((—k + 1)h)? 


) 


uf! 


+(k?(AC) i411 + (& = 1)?(AC)aga2 +++ + (B= (bk - 2))?(AC) a1 4 


Equations (2) and (4) give us 275-1 (AC) e443 + (AC)n41,.441 = 0. There- 
fore, the first term of the above truncation error, that is, the term with 
coefficient u;, is removed. We can obtain classes of the method, namely 
several orders of convergence, by utilizing the above truncation error and 
eliminating the coefficients of the various powers of h for different choices 
of a,’s and (;’s. However, since our goal is to obtain the highest order of 
convergence, so we choose aj’s and {;’s that have the following conditions: 


(i) They satisfy the following relation 
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k k 
Yik- (5 —1))?(AC)etag + 25 °(B)nsij + (B)ksipri = 9, 


j=l j=l 


which eliminates the second term of the truncation error, that is, the term 
with coefficient hu’. The above relation can be written as follows: 


= k 
ao +25) a; = Bo t+2)— B;, 
j=l 


j=1 


because from the definition of matrices A and B, we have: 


k k-1 
So(k- (7 -1)) 2(AC)k41,5 =-a9 -25 a, 
j=l j=l 

k k 
2) (B)asag + (Basser = Bo +27 Bj. 
j=l j=l 


In accordance with the papers related to splines, the following relation is 
provided for ag and {o: 


k-1 k 
a9 +25 >a; =1=fo+2)— 8; 
j=l j=l 
In other words, 
k-1 
a =1-25- aj, (6) 
j=l 
and 
k 
Bo =1-25- By. (7) 
j=l 


For more details, the interested readers are advised to see [11, 16, 20, 19, 
21, 22, 25) and other related papers. 


(ii) The remaining unknown elements, namely the following (2k — 1) co- 
efficients are chosen in such a way that the terms with coefficient hiu) 
pokey (t*) 

a 


to 
, in the truncation error, are eliminated: 


Q1,02,.--,Ak—1,F1, B2,.-., Br. 


Therefore the local truncation error associated with (5) is O(h***?), 
k = 1,2,.... Consequently, the proposed method is convergent of order 
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O(h**), k=1,2,.... 


2.3 Spline solution 


Now we discretize (1) as uy’ = fi+giui, 1 =1,2,...,n—1, at the grid points, 
where f; = f(x), 9; = g(a). Asa result, the vector V” can be rewritten as 


J 2k+1<j<n-1. 


o] 


oak n+ ats CSS eR +1, 


If we substitute V” in (5) for i=k,k+1,...,n—k, then a system with 
(n—2k+1) linear algebraic equations and (n—1) unknowns as uj, U2,..-, Un—1 
is obtained. Note that up = A and u, = ¥. It can be represented in the matrix 
form as follows: 


DU = R, (8) 


where U = [w1,u2,-.-,;Un—1]7 is a column vector with (n — 1) elements. 


Moreover, D is a matrix of order (n — 2k +1) x (n— 1), indicated as follows: 
D = (AC) +h? BG, 


with G = diag(g1,92,---;9n—1). Matrices AC and B are also defined by 
omitting the (k — 1) first and last rows of matrices AC and B, respectively. 
In other words, for 1 <i <<n—2k+1and1<7<n-—1, we have 


(AC)i,5 = (AC)e+u—y> (Boag = (B)eea—1,3- 


Finally, the vector Ris given by 


2k+1 


R= —W SO (Byers fj4i-2 
j=l 
—uo ((AC)e+1,1 +A? 90(B)e+11), = 1, 
a) a 2<i<n- 2k, 


—Un ((AC) 41,1 + h? gn (B)k-+1,1) ,ti=n—-2k+1. 


2.4 Development of the boundary formulas 


To obtain a unique solution for the system (8), we need (2(k — 1)) more 
equations; thus, we define them in the following form: 
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k+itl 4k-1 

y Glu; +h? » Bi uy = 0, i=1,2,...,k—1, 
j=0 

k+itl 4k-1 


dy Sn jun—j +h? 2 Pas ti =0, t=n—(k—1),...,n-—2,n—1. 


In order to use the band matrices in the new system, that is, system 
(8) along with the above equations, we use the following replacements for 
j=1,2,...,0+&: 

a, =(AC)ij, 1=1,2,...,k-1, 
ai =(AC)in-7, t=n—(k-1),...,.n-2,n—-1. 


These replacements simplify the pon vereerce analysis of the method. The 
other unknown coefficients, Bi ’s and a s, are determined by considering 
the local truncation error of order O(nAk+2) for the added equations and 
using Taylor’s expansion of these equations for i = 1,2,...,k —1 around 29 
(or for i= n—(k—-1),...,n—2,n—1 around z,). 


On the other hand, from (2 ”) and x ), we have (AC); 5 = (AC)n-in—j- 
Consequently, a = = ar- = and a= = ee i Considering these justifications, 
system (8) is onvanted to the following system: 


DEF, (9) 
with 
D=AC +h? BG, (10) 
where the matrix B is (n — 1) x (n — 1)-dimensional as the following form: 
(B)nin-j =Bi, 1<i<h-1, 1597 <4k-1 


(B)ig = (11) 
(B)i; for other i, j. 


For the column vector R with (n — 1) elements, we have 


uo (Gi, + h? go) — h? (Bi. fo + O57" (Big fia), 1<i<k-l, 
Ri = Ri-(e-1); k<i<n-k, 
un (ai, + h? gn fi,) — h? (Bi fn + 387 (B)igfn—j)) M-R+1Si<n-1. 


Finally, by solving the system (9), we obtain the solution vector U, the 
elements of which are approximately equal to the solution of (1) at nodes 
U1,U2,.++,en—-1- 
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3 Numerical results 


In order to test the viability of the proposed method and to demonstrate 
its convergence computationally, some BVPs including the cases of linear, 
nonlinear, perturbed, and system are considered. We measure the accuracy 
in the discrete maximum norm 


||| = ||U =. Vexact|l = max |U; = (Uexact)il; 


1<i<n-1 
and the convergence rate for linear and perturbed cases 


|B" | 
|r|” 


CR = log,( 


where ||E”|| and ||E?"|| are the maximum absolute errors on n and 2n grid 
points, respectively. The results are listed in tables for different choices of n 
and k. From the tables, we see that the quantity C'R is close to 4k for each 
k. In other words, by reducing the step size from h to A the observed errors 
are approximately reduced by a factor (Cas verifying the convergence order 
of the presented method, that is, O(h*"), k = 1,2,.... For example, in the 
rows related to n = 16 and n = 64 from Table (1), it is observed that the 
maximum absolute error is decreased by a factor (+)** when n = 16 is varied 
to n = 64. Namely, we have 


1 
fork=1: (6.72% (107°)) * Gy ~ 2.63 * (10~1"), 


1 
fork=2: (3.12% (1071°)) x Gr ~ 8.09 * (10~?"), 


The outcomes indicate that our presented method produces more accurate 
results in comparison with those obtained by other methods. It should be 
mentioned that the computations associated with the examples in this paper 
were performed using Mathematica 8.0. Run applications were done in just 
a few minutes. The numerical results in tables were written just for some 
values of k, but we could solve the examples for other values and the results 
are quite satisfactory as was already expected. 

In addition, we have used some plots to illustrate the behavior of the 
numerical solutions. Furthermore, since B is (n — 1) x (n — 1)-dimensional, 
in (11) we should have 4k — 1 < n—1 and k—1< n-—1, which results in 
4k <n. This can be seen in the results tables. 


Example 1. We consider the following linear two-point BVP: 


u" (x) — ula) = «7 — 2, 
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The exact solution is 


The corresponding maximum absolute errors and convergence rates in our 
computed solutions are listed in Tables (1) and (2), respectively. Rashidinia, 
Jalilian, and Mohammadi [21] solved this problem by using the nonpolyno- 
mial quintic spline. Although their method is similar to ours for k = 2, the 
only difference is that they used a lower order of convergence of the method, 
see Table (3). 

The graph of the exact and approximate solutions of Example (1) for 
n = 20 and k = 1,2,3,4,5 is depicted in Figure (1). 


| ! ! ! | ! ! ! | ! ! ! | ! ! ! | 
0.2 0.4 0.6 0.8 1.0 


Figure 1: Plot of the exact and numerical solutions of Example (1) for k = 1,2,3,4,5 
and n = 20 


Example 2. We consider the following nonlinear BVP, the classical Bratu’s 
problem: 


where 7 > 0. The exact solution is 
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where 0 = \/27cosh($). The Bratu’s problem has zero, one, or two solutions 
when 7 > ec, 1 = Ne, and 7 < me, respectively, where the critical value 7. 
satisfies the equation 1 = 4,/27.sinh(4) and it was evaluated in [5, 12] that 
the critical value 7, is given by 7, = 3.513830719. 

We have solved this example for 7 = 1,2, and 3.51 using our method with 
different values of & and tabulated the results in Tables (4), (5), and (6). Note 
that, in this example, we have used the Newton—Raphson algorithm, just with 
two iterations. Thus, there are errors related to the initial conjecture and 
the number of iterations, in addition to the error of our method. Tables (7) 
and (8) contain the comparison of our results and the results in [6, 13, 31]. 
The method in [31] is the same as our method for k = 2 with a lower order of 
convergence. Note that, the mentioned references have presented the outputs 
of their methods only for n = 10, so for the sake of comparison, in Table (7), 
we have to show the results of our method only for this value of n. We can 
use both k = 1 and k = 2 (according to condition 4k < n ), but k = 2 
provides better results. Thus, we display its maximum absolute error. 

Figure (2) plots the graphs of analytic and approximate solutions of Ex- 
ample (2) for n = 20, 7 = 1, and k = 1,2,3,4,5. 


0.12 


0.14 


0.10} 
0.08 } 
0.06 } 
0.04} 


0.02 F 


f f f f f dt f a a f 
0.2 0.4 0.6 0.8 1.0 


Figure 2: Plot of the exact and numerical solutions of Example (2) for k = 1,2,3,4,5 
and n = 20,7 =1 
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Example 3. We consider the following singularly Perturbed BVP: 


eu (x) = u(x) + cos?(rax) + 2en* cos(272), 
u(0) = u(1) = 0. 


The exact solution is given by 


—(1-2) —2x 
exp(—7-—) + exp(Fz) 
u(a) = We = ue cos? (Ta). 
1+ exp(=.) 


The maximum absolute errors and convergence rates for « = ie are tabu- 
lated in Tables (9) and (10), respectively. The results for this example from 
[4, 8, 17, 22, 27] are listed in Table (11). Note that the method used in [4, 22] 
is the same as our method for k = 2, but with a lower order of convergence. 


The results of [8] are also the same as the results of our method for k = 4. 


We observe from Figure (3) that the graphic of the approximate solution 
of Example (3) for n = 20 and k = 1,2,3,4,5 coincides with the graphic of 
the exact solution. 


Figure 3: Plot of the exact and numerical solutions of Example (3) for k = 1,2,3,4,5 
and n = 20 


Example 4. We consider a BVP in calculus of variations, that is, the prob- 
lem of finding the extremal of the functional [23]: 
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™ 


faa | © (ul? (a) + w2,(w) + 2uy(x)urr(a)) de, 


with boundary conditions 


eam ur(Z) =1 


urr(0) = 0, ur(S =-l. 


The exact solution is given by uy(x) = —uy7(x) = sin(x). For this prob- 
lem, the corresponding Euler-Lagrange equations are 


that is, a system of equations such as (1). It should be mentioned that in this 
example, we compute ||£,,,|| and || E.,,,|], but one of them is displayed in Table 
(12), because, the value of both norms is the same. This example has already 
been solved by using cubic [29] and quintic [30] parametric spline methods, 
namely, the same as our method for k = 1 and k& = 2 (with a lower conver- 
gence order), respectively. The sinc-Galerkin method [28] is also the other 
method that has been used for solving the above problem. The mentioned 
references have provided the numerical results only for n = 5,10, 20,..., 50. 
To make a proper comparison with these methods, we have shown our results 
only for these values, in Table (13). We have selected k’s that apply to con- 
dition 4k < n and give the best outputs. For instance, for n = 50, we could 
display the numerical results of our method for k = 1,2,3,...,12, but since 
k; = 12 gives the best result, we display its maximum absolute error. 

The numerical results of Example (4) for n = 20 and k = 1,2,3,4,5 are 
plotted in Figure (4). Note that we have displayed this graph just for u;(a). 
Similarly, it can be shown for uy;(z). 


4 Conclusion 


A long process is needed to obtain the differential relations of spline-based 
methods. Therefore, it is important to use a method that has considerably 
less computational effort with high accuracy and improves the spline meth- 
ods. In this paper, for the first time, a generalized form of methods based on 
parametric splines of degree (2k+1), k = 1,2,..., was introduced that has all 
the mentioned properties. A very good accuracy of this method was demon- 
strated for solving some linear, nonlinear, perturbed, and system of BVPs. 
We mention some advantages of our method in the following remarks. 


Remark 1. It is necessary to obtain the criterion of spline function in the 
spline methods. For instance, in the parametric spline method, this criterion 
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Figure 4: Plot of the exact and numerical solutions of Example (4) for k = 1,2,3,4,5 


and n = 20 
Table 1: Maximum absolute errors for Example (1). 
n k=1 k=2 k=3 k=4 k=5 k=6 k=7 
4 1.66E—6 - - 
8 107E—7 = 2.24E—12 = = - = - 
12 2.13E—-8 4.92E—-14 1.13E—18 = = - - 
16 6.72E—-9 3.12E—-15 2.48E—20 2.10E—25 = _ = 
20 2.76E-9 3.62E—-16 1.19E—21 4.44E—-27 1.74E—32 = = 
24 1.33E—-9 6.17E-17 9.97E—23 1.84E—28 3.64E—34 7.42E—40 - 
28 7.18E—-10 1.37E-17 1.21E—23 1.23E—29 1.34E—-35 1.52E—41  1.77E—47 
64 2.63E—11 8.09E—21 1.34E—28 5.35E—-36 2.28E—43 1.02E—50 4.75E—58 
128 1.64E—-12 2.61E—23  8.75E—-33 2.22E—41 6.07E—50 1.74E—58 5.16E—67 
256 = =1.02E—-13) 9.71E—26 5.52E—37 8.84E—47 1.52E—56 2.76E—66 5.18E—76 
512 6.48E—-15 3.74E—-28 4.03E—41 3.44E—-52 3.73E—-63 4.25E—74 5.01E—85 
1024 4.02E—-16 1.45E—-30 6.27E—45 1.32E—57 9.02E—70 6.43E—82 4.75E—94 
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Table 2: Convergence Rates, Example (1). 


n kK=1 k=2 k=3 k=4 k=5... 
64 4.00 8.27 13.90 17.87 21.84 
128 ~=4.00 8.07 13.95 17.93 21.92 
256 = 3.98 8.02 13.74 17.97 22.10 
512 3.99 8.01 12.65 17.99 21.83 


Table 3: Maximum absolute errors in [21] for Example (1). 


Second-order [21] 


Fourth-order[21] 


Sixth-order[21] 


8 1.09E—4 
16 3.06E—5 
32 8.11E—6 
64 2.09E—6 


5.22E—8 
2.31E—9 
1.34E—10 
8.42K—12 


8.75E—11 
5.74E—13 
2.30K—14 
3.68E—14 


Table 4: Maximum absolute errors for Example (2), 7 = 1. 


n ak k=2 k=3 k=4 ee 
4 1.17E—5 = 
8 7.23E—7 1.78E—9 a _ = 
12) 1.42E—-7 1.72E—-11 8.37E—-13 — _ 
16 4.50E—-8 5.50E—-13 6.32E—-15 5.89E—16 a 
20 1.84E-8 4.20E—-14 2.94F—-16 3.08E—-16 3.98E—16 
24 8.89E-9 6.63E—-15 2.22E—-16 2.77E—-16 2.91E—16 
28 4.79E-9 1.38E—-15 4.16E—-16 4.99EF—-16 2.49E—16 
32 2.81E-9 1.41E—-15 5.82E—-16 4.30E—-16 4.44E—16 
36 «61.75E-9 3.33E—-16 1.05E-15 6.10E—-16 1.38EK—16 
Table 5: Maximum absolute errors for Example (2), 7 = 2. 
n ek k=2 k=3 k=4 ee 
4 1.58E—4 a 
8 9.55E-6 1.22E—7 _ — a 
12 1.87E-6 4.09E—-10 3.32E—10 _ _ 
16 5.92E—7 8.14E—-11 2.38E—-12 1.07E—12 _ 
20 2.42E—7 1.09E—-11 3.06E—-13 1.78E—-14 3.88E—15 


Iran. J. Numer. Anal. Optim., Vol. 13, No. 4, 2023,pp 578-603 


Sarvari 598 
Table 6: Maximum absolute errors for Example (2), 1 = 3.51. 
n k= 1 k=2 k=3 k=4 k=5 
4 5.505K—1 = 
8 2.88E—3 2.29E—4 = - 
12) 5.51E-4 2.46E—-5 9.43E—6 a _ 
16 1.73E—-4 6.92E—7 9.35E—-7 4.08E—7 a 
20 7.08E—5 1.18E—-8 4.85E-9 4.00E—-8 1.74E—8 
Table 7: Comparison of ||E|| for Example (2), 7=1 n= 10. 
nm Our method fork =2 Method[31] Method|{6] 
10 1.44E—10 5.87E—10 8.89E—6 
Table 8: Comparison of ||£|| for Example (2) with 7 = 1,2 and 3.51. 
Our method Method[13] 
n n=l n=2 n = 3.51 n=l n=2 n = 3.51 
8 178b—9k—2) 122E—7(k=2)  2.29b—4(k— 2) 5.04b—-9 453-8  3.51B—5 
16 5.89E-16(k=4) 1.07E—12(k=4) 4.08E—7(k= 4) 4.66E-11  1.76E-9 1.45E—7 
32 4.30E-16(k=4) 6.71E-15(k=3) 9.32E-10(k=2)  8.33E-13  2.13E-11 —‘1.02B-9 
64 471E-16(k=6) 3.60E—-15(k=3)  1.37E—9(k= 2) 9.21E—15 2.87E-13  1.48E—11 
128 2.292B-15(k=6) 4.99B-16(k=6) 1.37E—9(k= 2) = 247B-14  1.58E-13 
Table 9: Maximum absolute errors for Example (3), ¢ = +. 
n k= k=2 k=3 k=4 k=5 k=6 k=7 
4 1.15E—2 a a 
8 6.65E—4 4.45E—5 = = = = = 
12 1.29E-4 7.49E-8  5.02E—8 - = - a 
16 4.07E—5 2.73E-8 4.75E—-10  1.72E—11 - - - 
20 1.66E—5 5.10E—-9 1.50E-12 2.17E-13 2.45E—15 - = 
24 8.01E—6 1.05E-9 1.03E—-12 3.93E—-15 3.53E—-17 1.73E—19 a 
28 4.32E-6 2.58E-10 2.17E-13 1.51E-17 7.65E—-19 2.67E—21 6.76E—24 
64 1.58E—-7 9.10E—-14 4.72E—-18 2.66E-—22 1.51E—26 8.18E—31 3.78E—35 
128 9.87E-9 2.26E—-16 3.28E—22 1.28E—27 5.36E—33 2.31E—38 1.02E—43 
256 6.17E—-10 9.47E—-19 2.08E—26 5.21E—33 1.40E—-39 3.95E—46 1.14E—52 
512 3.86E—11 3.76E—21 1.29E—-30 2.02E—38 3.43E—46 6.10E—54 1.12E—61 
1024 2.41E—-12 1.47E—-23 7.92E—35 7.78E—44 8.25E—53 9.18E—-62 1.05E—70 
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Table 10: Convergence Rates, Example (3). 
n k=1 k=2 k=3 k=4 k=5... 
64 4.00 8.65 13.81 17.66 21.42 
128 =3.99 7.89 13.94 17.90 21.86 
256 3.99 7.97 13.97 17.97 21.96 
512 ~=4.00 7.99 13.99 17.98 21.98 
Table 11: Comparison of |||] for Example (3), ¢ = ‘re 
n Method|8] Method{17| Method[4] Method{22] Method[27] 
16 1.72E—11 1.22E—6 1.57E—5 4.07E—5 1.20EK—4 
32 1.52E—17 6.45E—9 8.79E—7 2.53E—6 7.47E—6 
64 2.66E—22 3.40E—11 5.32E—8 1.58E—7 4.67E—7 
128 = 1.28E—27 1.03E—12 3.30E—9 9.87E—9 2.90E—8 
Table 12: Maximum absolute errors for Example (4). 
n k=1 k=2 k=3 k=4 k=5 k=6 ae 
4 2.76E—5 = = - a - - 
8 1.72E—6  1.80E—10 - - - - - 
12 3.41E—7 3.37E-12 5.75E—16 - - - - 
16 1.08E—7 1.97E—-13 1.03E—-17 6.79E—22 = = = 
20 4.44E—-8 2.32E-14 4.72E—-19 1.22E—23 3.53E—28 aa = 
24 2.14E—-8 4.16E-15 3.74E—-20 4.60E—25 6.41E-30 9.33E—35 = 
28 1.15E-8 1.00E—-15 4.37E—21 2.92E—26 2.16E—31 1.69E-36 1.37E—41 
512 1.03E-13 3.66B—26 1.97E-38 5.96E—49 3.95E—59 2.75E-69 1.98E—79 
Table 13: Comparison of ||E|| for Example (4). 
n Our method Method[30]_ Method[29} Method[28} 
5 1.12E—5(k= 1) — 1.12E—5 — 
10 = 2.02E—11(k= 2) 6.70E—10 7.05K—7 2.72E—4 
20 =3.53E—28(k= 5) 7.07E—12 4.44 —8 8.69E—6 
30»: 1.74E—42(k= 7) 8.10E—13 8.77E—9 5.65E—7 
40 2.61E—63(k= 10) 1.55E—13 2.78E—9 5.47E—8 
50 9.59E—67(k= 12) 4.21E—14 6.93E—9 
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is obtained by solving a special ordinary differential equation. Or nonpolyno- 
mial spline is a function with unknown coefficients that should be determined 
accordingly. In all these cases, some time-consuming calculations are needed, 
while the criterion and coefficients of no function are required in our method. 


Remark 2. The continuity property of spline and its derivatives in grid 
points plays a major role in all of the spline methods. One can use this 
property to obtain the required spline relations. In this paper, instead of 
using the properties of spline directly, to save time and reduce calculations, we 
derived the consistency relations from a special algorithm and then obtained 
its matrix form by defining some band matrices. 


Remark 3. The approximate solution converges to the exact solution of 
order O(h4*). It follows that ||E|| + 0 as h + 0. The convergence occurs 
more quickly when k is a larger number. Indeed, the order of error is not 
fixed and decreases by increasing the value of k. It is regarded as one of our 
method’s advantages. In addition, since we have h = bea it concludes that 
|Z || = O((=*)**). This indicates that an increasing k is more effective than 
that n in reducing error. It can be seen in the tables containing numerical 
results. 


Remark 4. We claim that the proposed method can be applied to solve other 
similar differential equations in particular as u?™ (ax) = f(x) + g(x)u(2), 
where m indicates a positive integer. This will be considered in our future 
research. Moreover, because of the adequate flexibility and expandability of 
this method, there is a possibility of achieving the generalized form of the 
methods based on splines with the degree (2k), k =1,2,.... 
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